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Abstract—In this work, we show how we can improve 
the image resolution capabilities of a Phase Conjugating 
(PC) lens as well as the angular resolution of Luneburg 
lens antennas, aperture antennas, and the Randomly 
Distributed Radar Array (RDRA), by employing signal 
processing techniques, such as the Correlation Method 
(CM), the Minimum Residual Power Search Method 
(MRPSM), the sparse reconstruction method, and the 
Singular-Value-Decomposition (SVD)-based basis matrix 
method. In the first part, we apply these techniques for 
sub-wavelength imaging in the microwave regime by 
combining them with the well-known phase conjugation 
principle. We begin by considering a one-dimensional 
microwave sub-wavelength imaging problem handled by 
using three signal processing methods, and then we move 
on to two- or three-dimensional problems by using the 
SVD-based basis matrix method. Numerical simulation 
results show that we can enhance the resolution 
significantly by using these methods, even if the 
measurement plane is not located in the very near-field 
region of the source. We describe these proposed 
algorithms in detail and study their abilities to resolve at 
the sub-wavelength level. Next, we investigate the sparse 
reconstruction method for a normal Luneburg lens 
antenna and the Correlation Method and the SVD-based 
basis matrix method for a flat-base Luneburg lens 
antenna to estimate the Direction-of-Arrival (DOA). 
Numerical simulation results show that the signal 
processing techniques are capable of enhancing the 
angular resolution of the Luneburg lens antenna, enabling 
the lens to locate multiple targets with different scattering 
cross-sections, and achieving higher angular resolution. 
Finally, we present a hybrid approach which combines 


the Correlation Method and the SVD-based basis matrix 
method to achieve considerably higher angular resolution 
when using two representative aperture antennas, e.g. a 
parabolic reflector, a random array, to estimate the DOA. 
Numerical results show that the proposed hybrid 
approach can help achieve a higher angular resolution 
than that without the use of signal processing. 


Index Terms—Phase Conjugating (PC), Direction of 
Arrival (DOA), Randomly Distributed Radar Array 
(RDRA), Correlation Method (CM), Minimum Residual 
Power Search Method (MRPSM), Singular Value 
Decomposition (SVD). 


I. INTRODUCTION 


The problem of achieving high resolution focusing by using 
a “superlens” or a “perfect lens” [1]-[2] has recently attracted 
considerable attention of researchers, since such lenses have 
the potential to overcome the diffraction limit and achieve 
superior resolution in microwave imaging and in angle 
estimation type of problems. 


In microwave imaging, the topics of Phase Conjugation 
(PC), as well as Time Reversal Mirrors (TRMs) have been 
thoroughly investigated and a great deal of research has been 
conducted [3]-[4]. Microwave imaging with a resolution at 
the sub-wavelength level is attractive in many applications, 
which include: high-directivity arrays used in modern 


self-tracking wireless communication systems [5]; 
electromagnetic sounding of 3D_ structures of an 
inhomogeneous dielectric half-space for nanophysics, 


biological and medical diagnostics [6]; and short-range 
Ultra-wide Band (UWB) radar imaging [7], to name a few. 
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Many approaches have been proposed for microwave 
sub-wavelength imaging, including nonlinear metamaterial 
elements [8], periodic layered metal-dielectric structures [9], 
metallic screens comprised of closely spaced and unequal slits 
[10], Fresnel zone-plates [11], etc. Several devices based on 
PC lens have been reported, e.g., a phase conjugating lens 
consisting of a double-sided assembly of straight wire 
elements [12], metallic strip gratings which perform 
evanescent-to-propagating wave conversion [13], 
sub-wavelength array of planar monopoles [14] and split-ring 
resonators loaded with varactor diodes [15], etc. 


However, the applications of the current techniques still 
pose huge challenges. For instance, to realizing a “superlens” 
by using artificial materials or matematerials not only poses 
many practical difficulties, but complex is costly as well. 
Moreover, a common feature of these existing imaging 
techniques is that they must rely on the availability of 
near-field information, which is not too convenient to acquire 
in practice. Recently, signal processing techniques have been 
proposed into lens imaging to achieve improved resolution 
[16]-[17]. In [18]-[22], signal processing techniques for 
sub-wavelength imaging in the microwave regime by 
combining them with the PC lens, namely the Correlation 
Method (CM) [18], the Minimum Residual Power Search 
Method (MRPSM) [19]-[20], the sparse reconstruction 
method [18], and the Singular-Value-Decomposition 
(SVD)-based basis matrix method [21]-[22]. We demonstrate 
in this work that these techniques can achieve sub-wavelength 
resolution even when the measurement plane is not located in 
the near-field region of the source. The proposed techniques 
are detailed below and their abilities to resolve at the 
sub-wavelength level are examined. 


In the Direction-of-Arrival (DOA) estimation problem, 
wide-angle scanning and high-angular resolution are often 
desired in many applications, e.g., wireless communication, 
air navigation, automotive radar, etc., [23]-[24], and 
mechanically steerable antennas are traditionally used for 
such applications. However, such antennas are usually 
operated by a mechanical rotation type of device which 
cannot meet the requirements of rapid beam-scanning and 
flexible control. It is well-known that phased-arrays are more 
sophisticated and can serve multiple functions for different 
applications. However, they can be bulky, expensive to 
fabricate and often require signal-processing hardware. The 
Luneburg lens [25], which has the capability of all-angle scan 
regardless of the frequency of operation, as well as excellent 
focusing characteristics, is an attractive candidate for many 
applications such as multi-beam antennas, multi-frequency 
scanning, and spatial scanning [26]-[27]. Many researchers 
have focused on the problem of manufacturing and realization 
of the Luneburg lens, including the design and 
implementation of 2D or 3D Luneburg lens antennas. 
Unfortunately, its applications in the area of DOA estimation 
are rarely mentioned. 


In [28], Xin et al. reported an approach for detecting a 
single target by using five conformal detectors located on the 
curved surface of a single Luneburg lens antenna to receive 
the signal from -20° to 20°. The operating frequency is 5.6 
GHz, and the initial direction finding results from using a 


correlation algorithm showing that the estimated error is 
smaller than 2° within the angular range of incident angles 
from -15° to 15°. However, the range of incident angles is still 
limited, and the number of detectors (5) is not adequate to 
achieve the desired level of accuracy. Recently, we have 
proposed a novel flat-base Luneburg lens antenna, and have 
investigated its scan performance and compared it with 
conventional lenses [29]. Based on the previous studies, we 
propose the use of a Luneburg lens antenna for DOA 
estimation by using signal processing techniques, i.e., the 
sparse reconstruction method, the CM [30], the SVD-based 
basis matrix method. Numerical results show that, we can not 
only achieve both wide-angle scanning and high-angular 
resolution by utilizing signal processing techniques, we can 
also operate over a broad frequency range and for different 
polarizations. 


Estimating the DOA by using Uniform Circular Arrays 
(UCAs) or Uniform Linear Arrays (ULAs) is widely used in 
communication and radar fields [31]-[32], and a number of 
angular resolution algorithms, e.g., Multiple Signal 
Classification (MUSIC) [33], and Estimation of Signal 
Parameters via Rotation Invariance Techniques (ESPRIT) 
[34], have been developed for different applications. However, 
these techniques do not perform too well when applied to 
highly correlated or coherent signals generated via multipath 
propagation scenario. Moreover, the DOA estimation carried 
out by using the UCAs or the ULAs require strict geometry 
modeling and amplitude and phase calibration of the sensors; 
furthermore, they become unreliable when some of the 
sensors fail, and the array is no longer periodic. To circumvent 
these problems, a hybrid approach which combines the CM 
and the SVD-based basis matrix method for DOA estimation 
for aperture antennas, for instance a parabolic reflector, has 
been proposed to improve the estimation in a computationally 
efficient manner. The proposed method uses the correlation 
approach to estimate the DOA within a certain angular range. 
The inverse matrix method and the SVD-based basis matrix 
method are subsequently employed to improve the resolution 
of the aperture antenna. Furthermore, we investigate a 
framework based on the Randomly Distributed Radar Array 
(RDRA) for DOA estimation, and also present a hybrid 
approach that combines the correlation method and the basis 
matrix method with SVD to achieve substantially improved 
angular resolution over that is achieved by either of the two 
methods if applied individually. The hybrid method involves 
generating a data matrix using measurements performed for 
different incident angles by recording the outputs of each 
sensor of the array, obtained either from simulations or real 
measurements. Following this, we correlate the measured data 
with the generated data matrix, to derive an initial estimate of 
the angular ranges of probable DOAs. We then apply the basis 
matrix method based on the SVD of the data matrix to 
estimate the angles of incidence accurately. Numerical results 
show that the proposed hybrid approach can achieve a high 
resolution even when the sensors are randomly distributed in 
the array. 
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II. MICROWAVE SUB-WAVELENGTH IMAGING WITH A 
PHASE CONJUGATION LENS 


A. One-Dimensional Microwave Sub-wavelength Imaging 
(1) Principle of Phase Conjugation Lens 
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Fig. 1. Scheme for microwave sub-wavelength imaging via PC. 


A representative scheme of microwave sub-wavelength 
imaging via the use of Phase Conjugating (PC) lens is shown 
in Fig. 1. The test objects are y-oriented Perfect Electric 
Conductor (PEC) strips illuminated by a plane wave, that are 
located in the x-y plane (source plane) at z = 0. The scattered 
field distribution is first measured on a grid set up on the 
measurement plane at z = d. Next, the measured data is 
phase-conjugated and propagated to the image plane at z = 2d 
to obtain an image of the original object in the image plane. 


The spatial field distributions scattered by the object, 
sampled at the measurement and image planes, can be related 
by using a transform method. The Fourier transform of the 
spatial field distribution over the x-y plane located at z = 0 can 
be expressed as: 


F(k,,k,,z=0)= | [raz = (ele! dxdy (1) 
where k, and k, denote the spectral variables; f(x, y, z = 0) and 


F(k,, ky, z = 0) represent the spatial and spectral distributions 
at the source plane, respectively. 


Next, we derive the spectrum G(k,, k,, z = d) at the 


measurement plane as follows by propagating F(k,, k,, z = 0) 
to a distance d along the z-direction: 
G(k,,k,,z=d)=F(k,,k,,z=O)e (2) 


The wavenumber k, in the z-direction is related to ką, k, and 
k (k = 27/2) as follows: 


he k? — (k° +k?) when k? +k? + à (a) (3) 


© J-j i(k? +k?)-k? when k2+k2>k? (b) 


The form of k,, as expressed by (3a), corresponds to the 
propagating waves (visible part of the spectrum), while (3b) 
defines the evanescent waves (invisible range of the 
spectrum). 


The next step involves phase conjugation followed by 
propagation of the spectrum G(k,, k,, z = d) to the image plane 


to derive H(k,, k,, z = 2d). This is accomplished by using the 
expression shown below: 


A(k,,k,,Z=2d) = G'(-k, -k,,2 = dye it 
= F’(-k,,-k,,z = 0) 


(4) 


which relates the spectrum of the distribution at the image 
plane to that at the source plane. 


After having performed the operations of Fourier transform, 
phase conjugation, and propagation to the image plane, we 
find that the spectrum at the image plane is a “flipped” (with 
respect to the spectral variables k, and k,) and conjugated 
version of the original spectrum at the source plane. By 
applying the well-known properties of Fourier transforms, the 
relationship between the spatial field distribution at the source 
and image planes can be expressed as follows: 


h(x, y,z=2d) = f (x,y,z =0) (5) 


It is evident that the spatial field distribution at the image 
plane is the conjugate of the original spatial field distribution 
at the source plane. 

In general, the resolution of the lens is always 
diffraction-limited if the phase conjugating lens is positioned 
in the far-field region of the sources, because the evanescent 
component of the spectrum carrying sub-wavelength 
information cannot reach the phase-conjugating surface. 


(2) Correlation Method 


Next, we consider the problem of generating an image for 
an arbitrary combination of PEC strips, whose number is 
unknown as yet, and whose interelement space is in the 
sub-wavelength regime. 


Correlation coefficient 


1 source 


2 sources 


[1] = 


a Peak value = | | | 


4 sources 


N sources 


Different combinations 
Fig. 2. Representative scheme of CM. 


We now begin to describe the Correlation Method (CM), 
which is widely used in signal processing. As shown in Fig. 2, 
we assume that an arbitrary number of sources are present at 
the source plane, ranging from a single source, to a maximum 
number N. For each combination, we calculate the radiated 
fields at the measurement plane and then correlate them with 
the measured data. The peak value of the correlation 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


coefficients is used to determine the correct position of the 
original sources at the source plane. 

The CM is based on correlating the field distribution 
generated by each of the assumed combinations of sources 
with the original one. Without loss of generality, we assume 
that we are dealing with the £\-component of the electric field. 
The correlation coefficient is defined as: 

1 
(MM ; -louc, 


M; M ; 
YSE" any) -Etan y- ua) | 


i=l j=l 


Corr.coef(n) = 


(6) 


where Uy and wu, represent the mean value of the field 
distribution of the assumed combinations of sources 
(subscript 'A') and the measured field (M), while oy and o4 
are their standard deviations; respectively. M; and M; are the 
numbers of samples along x and y gathered at the 
measurement plane. Here the correlation coefficients are 
complex; hence, we compare their magnitudes. 


Specifically, all of the combinations are generated in a 
search area first identified by using the image obtained via an 
application of the phase conjugation algorithm. The search is 
not carried out at locations where the values of the normalized 
image fall below a certain threshold level, say 0.5. 
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Fig. 3. Correlation coefficients (left column) as a function of the 
considered combination of sources when the separation distance is (a) 24/5, 
(c) A/5 and (e) A/10 ; and normalized spatial images (right column) obtained 
by PC (solid line) and CM (circle marker) when the separation distance is (b) 
2A/5, (d) A/5 and (f) 2/10. 


The CM can achieve a very high resolution if one can 
accurately locate the maximum correlation coefficient. 
However, the CM approach has some limitations of its own. 
First, the maximum correlation coefficient used to determine 
the correct source combination is only slightly larger than the 
one which is the next highest, especially when the sources are 
close to each other or the Signal-to-Noise Ratio (SNR) is not 
high. Secondly, during the process of searching for the 
sources, the number of combinations to be considered can 
become very large and the computational burden can increase 


4 


dramatically. For example, when the sources are randomly 
distributed in 10 positions, the number of combinations is 2 
however, it rises to 2° when the source distribution involves 
20 unknowns. 


In the discussion presented above, we have assumed that the 
magnitudes of the sources are all equal. If we wish to include 
sources With different strengths, we must consider an 
additional degree of freedom, which will further increase the 
number of required operations, and will negatively impact 
upon the computational efficiency of the CM. 


For the CM, we set the operating frequency to be 10 GHz, 
the spatial sampling rate at the measurement plane to be 1/20 
and the number of samples as 801 and the aperture size is 404. 
We assume that the distance d between the source and 
measurement planes as 2. We consider two sources centered at 
the origin along the x-direction, and assume that their 
separation distance 1s 24/5, A/5 and A/10; respectively. 


Figure 3 shows the variations of the correlation coefficients 
(left column) as a function of the various combinations of 
sources considered when the separation distances are (a) 24/5, 
(c) A/S and (e) 4/10. The spatial images (right column) are 
derived by using the phase conjugation and the CM when the 
separation distances are: (b) 22/5; (d) 2/5; and (f) 2/10. In the 
CM, we search for the maximum correlation coefficient, 
which corresponds to the correct combination of the sources. 
The present example shows that for a moderate number of 
sources, the CM performs well and is able to achieve 
sub-wavelength resolution, even with separation distances as 
small as 4/10. 


(3) Minimum Residual Power Search Method 


To circumvent some of the issues encountered in section 
above with the CM when the number of sources is not small, 
we apply the Minimum Residual Power Search Method 
(MRPSM), which systematically adjusts the local distribution 
of the sources by using an iterative minimum residual power 
search algorithm to obtain both the numbers and the positions 
of the sources. 


The step-by-step procedure for this method is as follows: 


1) Generate the image result from the measurement data 
by using the PC method. 


2) Extract the y 
normalize the curve. 


3) Determine the image region, N7p, by using the criterion 
where this region only covers the locations where the level of 
the normalized field is higher than a preset threshold ôV, e.g., 
OV =0.5. 


4) 


5) Find the position of the maximum value of the 
normalized y = O cut curve. 


6) Adjust the local source distribution in the 
neighbourhood of this position with the local adjustment 
region N;,4r, and calculate the image field of each source 
distribution. 


7) Search the source distribution which has the minimum 
power of the residual between the image field generated by 
the measured data and that calculated by the assumed source 


O cut from the image result and 


Set the local adjustment region to be N; yp. 
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distribution, according to the Minimum Residual Power 
Criterion (MRPC). 


8) Repeat until the search is finished and all the 
combinations of the source distribution within the image 
region Nr have been considered. 


9) Take the distribution which has the minimum power of 
the residual as the final result. 


The MRPC is defined as: 


amy 
A(n) = min SSE" ay- E] (7) 


i=l j=l 


where E” represents the field distribution of the measured 
field (M^, and E* represents the field distribution of the 
assumed combinations of sources (subscript 'A’) at the n-th 
search step, in which the sources are rearranged only in the 
local adjustment region. The quantity A(n) represents the 
minimum power of the residual at the n-th step. 
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Fig. 4. Schematic representations of (a) SMBFs and (b) SRBFs. 
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Fig. 5. (a) E, field distribution at the measurement plane; (b) image 
recovered by PC at the image plane; (c) normalized magnitude of E, at the 
image plane (y = 0 cut) with a preset threshold, i.e., 0.5; and (d) original 
source distribution; black strips represent the PEC strips. 


In both the CM and MRPSM approaches, we need to 
calculate the field distributions emanating from different 
combinations of these objects, and then correlate or compare 
them to the measured field. To numerically model the current 
distribution over the sources, we employ either the Sinusoidal 
Macro Basis Functions (SMBFs) [35], or the Sinusoidal 
Rooftop Basis Functions (SRBFs) [36], although we could 
equally use other basis functions as well. 


As shown in Fig. 4, SMBFs assume that the induced surface 
current on the PEC strips is sinusoidal over the entire length of 
the dipole, while when using SRBFs, the current distribution 
iS expressed as a superposition of these basis functions. 
Typically, when using SRBFs, a half-wavelength strip whose 
length and width are 1/2 and 1/20, respectively, is divided into 
10 and 3 sections in the x- and y-directions, respectively, and 
this makes the number of unknowns of the current distribution 
to be 27 times larger than that when we use the SMBFs 
method. We find that the SRBFs yield results are more 
accurate than those obtained with the SMBFs, because of the 
current distribution model is more accurate. 


Numerical simulations of different distributions of PEC 
strips have been carried out by using the FEKO 6.2 Suite [37] 
to generate the “measured” data. Due to physical constraints, 
the aperture size cannot be chosen as large as desired. In the 
following examples, we set the operating frequency to be 1 
GHz, the distance between the source and measurement 
planes to be 5A, and the observation aperture size to be 204. 
The equivalent sources are y-oriented PEC strips, each with 
length 4/2 and width 4/20, respectively. The sampling interval 
is A/S, both in the x- and y-directions, while the corresponding 
number of samples in the measurement plane is 101. 
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Fig. 6. Recovered results by (a) CM with SMBFs; (b) CM with SRBFs; (c) 
MRPSM with SMBFs; and (d) MRPSM with SRBFs; black strips represent 
the PEC strips. 


For the first case, 6 PEC strips are located at the following x 
positions: (-4, -3, 0, 1, 2, 4)x4/10, as shown in Fig. 5(d), and 
they are illuminated normally by a TE-polarized plane wave. 
Fig. 5(a) shows the electric field distribution observed over an 
aperture size of 204 x204 at the measurement plane, and it is 
simulated by using FEKO 6.2 Suite. Fig. 5(b) shows the 
image generated by utilizing PC only, and it is apparent that 
this technique alone is unable to resolve the position of the 
sources, even though the maximum interval between the 
sources is 34/10. Fig. 5(c) displays the y = O cut of the 
normalized field distribution at the image plane. We set the 
threshold ôV to be 0.5, and find that the region in the image 
plane where the values of the normalized field is above ôV = 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


0.5 ranges from -24/5 to 24/5. We set the unit interval to be 
A/10, so this image region contains 9 unknown units, implying 
that Nir = 9. In the CM, the number of combinations is 512. In 
the MRPSM, we set N,4r8 to be 9. 


Although the SMBFs method is 8 times faster, it can be seen 
from Figs. 6 (a) and (c) that the results are not very accurate 
when we use SMBFs in the CM or MRPSM. In contrast to this, 
it can be seen from Figs. 6(b) and (d) that the images are 
faithful reproductions of the object (array of strips) when we 
use SRBFs in the CM or MRPSM, albeit at an increased 
computational cost. The reason is that the SRBFs provide a 
better numerical representation of the current distribution 
over the strips. 


Next, we consider the case of 15 PEC strips that are located 
at the following x-positions: (-20, -15, -10, -9, -8, -7, -6, 0, 3, 
12, 13, 14, 15, 16, 20)xA/10, respectively, as shown in Fig. 
7(d). These 15 PEC strips are illuminated by a TE-polarized 
plane wave, whose incidence angle is (9 = 45°, ọ = 0°). The 
minimum separation distance between the strips is 4/10. The 
image region is larger than 4/, as may be seen from Fig. 7(c); 
hence, the CM would not be as efficient in this case as it was 
in the previous example, even if we use parallel computing 
techniques to help reduce the computational burden to some 
extent. However, due to its inherent characteristics, the 
MRPSM is suitable for the case where the image region 1s 
wide, even though its computational complexity increases as 
we increase the width of the image region. 


From the previous discussion, we see that the SMBFs have 
better efficiencies than the SRBFs, while the SRBFs have 
better accuracy in comparison to the SMBFs. In view of this, 
we propose a combination scheme, which begins with 
MRPSM with SMBEFs for the initial estimation and follows 
this up with the final refinement by utilizing a combination of 
the MRPSM with SRBFs. Thus, combing the advantages of 
both the SMBFs and SRBFs methods, this approach can not 
only generate the preliminary results faster, but can also refine 
it to achieve a more accurate and reliable result with greater 
efficiency. 


Figure 7(a) shows the results for the simulated E, field 
distribution derived by using FEKO, observed over an area of 
202 x204 in the measurement plane. Fig. 7(b) shows the 
recovered image generated directly by utilizing the PC, and it 
shows that the distribution of the sources is comprised of 
several clusters. Fig. 7(c) shows the y = O cut of the 
normalized field distribution at the image plane, and the 
image region ranges from -224/10 to 234/10 when we set the 
threshold ôV to be 0.5. This image region contains 46 
unknown units, 1.€., N; = 46; and we set Nzag to be 9. Fig. 7(e) 
shows that the MRPSM applied by using the joint method 
generates a final recovered result, which is a faithful 
reproduction of the original distribution, as shown in Fig. 7(d). 
In addition, we compare the running times of the MRPSM 
when using the SMBFs, SRBFs, as well as the joint method, 
and we find that the MRPSM with SMBFs is 6 times faster 
than that with SRBFs, and the MRPSM with joint method is 2 
times faster than that with SRBFs. In conclusion, the 
combination of MRPSM with SMBEFs is the fastest approach, 
but it does not produce an accurate result. In contrast to this, 


the SRBFs and the joint method can recover the image results 
with good fidelity, although at an increased cost. To 
summarize, the MRPSM _ offers good combination of 
efficiency and accuracy when combined with the joint 
method. 
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Fig. 7. (a) E, field distribution at the measurement plane; (b) image 
recovered by PC at the image plane; (c) normalized magnitude of E, at the 
image plane (y = 0 cut); (d) original source distribution; and recovered 
results by (e) MRPSM with SMBFs, (f) MRPSM with SRBFs, and (g) 
MRPSM with joint method; black strips represent the PEC strips. 
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Fig. 8. E, field distribution at the measurement plane when (a) SNRy = 20 
dB, (b) SNRy = 10 dB; recovered results by MRPSM with joint method when 
(a) SNRy = 20 dB, (b) SNRy = 10 dB; black strips represent the PEC strips. 


We now turn to the issue of the SNR of the measured data. 
Fig. 8 shows the numerical results for different levels of SNR 
of the measured data (SNRy) in the presence of additive 
Gaussian white noise. Figs. 5-8 demonstrate that both the CM 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


and the MRPSM can achieve a A/10 sub-wavelength 
resolution relatively efficiently when the image region is not 
wide. On the other hand, for the case when the image region is 
wide, the joint version of the MRPSM method, which 
combines the advantages of both SMBFs and SRBFs, 
produces accurate and reliable results with good efficiency. 
Specifically, such a joint approach can reliably achieve a 
sub-wavelength resolution of 4/10 even when the image 
region 1s wide. 

It is important to point out that in practice the local 
adjustment region N;,4r needs to be carefully chosen for a 
complex model. For instance, a larger value of N;4r May 
unnecessarily increase the computational complexity, while a 
value of N;4rR which is too small could reduce the accuracy of 
the results, because the mutual coupling effects may not be 
modeled sufficiently accurately. 


(4) Sparse Reconstruction Method 


The microwave sub-wavelength imaging problem can be 
treated as a sparse reconstruction problem when the 
distribution of the targets is sparse, and it may be stated as 
follows: find the positions and amplitudes of the excitation 
sources from the measured data and the characteristics of the 
given system. Here, we propose the some sparse 
reconstruction algorithms for this problem that are used in 
Compressive Sensing (CS) [38]. 


The construction of the microwave sub-wavelength 
imaging model is essentially based on the fields radiated by 
objects acting as secondary sources, e.g., different types of 
dipoles or perfectly conducting strips that are illuminated by a 
plane wave. Without loss of generality, we assume that the 
objects we are dealing with are half-wavelength dipoles. The 
induced current distribution on such dipoles can be 
approximated as follows: 


aahi] (eest 


I(y)= (8) 


| l 
al. sin} k} — + ——< y<0 
Fe ce | f | 


where J,, = constant is the magnitude of the exciting current, / 
is the length of the dipole. 


According to [39], we can calculate the £, field distribution 
as follows: 


— jkR, — jkR, — jkr 
B ea ee ole (9) 
á 4xz| k R, 2A T 


where R; and R, are the distances from the end points of the 
dipole to the observation point, and r represents the distance 
from the central point of the dipole to the observation point. 











We can establish the relationship between the measured 
data and the current distribution of the sources in a matrix 
form as: 

E=A.-I (10) 
where E represents the measured data vector which is 
consisted of E,-component, A is the system matrix and J is the 
source distribution vector. 
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Fig. 9. (a) Original distribution of strips; (b) image recovered by LSM; and 
(c) image recovered by SLO. 


We employ SRBFs to simulate the field distribution 
scattered by 10 PEC strips located along the x-direction, 
whose lengths and widths are 4/2 and //20, respectively. A 
plane wave, which propagates perpendicularly to the source 
plane with its electric field polarized along the y-direction, is 
incident upon these strips. We set the operating frequency to 
be 10 GHz. We assume that the induced current is sinusoidal 
in the y-direction and constant along the x-direction. The 
scattered electric field is measured at a distance d = À from the 
source plane, the spatial sampling rate at the measurement 
plane is 4/10, the number of samples is 101 and the aperture 
size is 104. The strips are centred at 
(-10, -8, -5, -4, -3, 0, 3, 5, 7, 9)x4/10 along the x-direction, 
respectively. 


Figure 9(a) shows the original distribution of the strips, 
where the bright bands represent the PEC strips with different 
induced current magnitudes. Instead of, Fig. 9(b) shows the 
image recovered by the Least Square Method (LSM), while 
Fig. 9(c) shows the recovered image obtained by using the 
Smoothed Jp algorithm (SLO) [40]. The least square method is 
a widely used approach for solving inverse problems. 
However, it is difficult to achieve an accurate solution by 
using this method, since the condition number of the system 
matrix A is very large, and the rank of A is much less than the 
number of unknowns. For this reason, we use the smoothed lo 
algorithm, which was originally developed for sparse 
reconstruction in CS, to find the sparse solutions of an 
underdetermined system of linear equations. 


Although the performance of this algorithm is good, it is 
sensitive to noise because of the ill-conditioned nature of the 
system matrix Á and its final optimization results depend on a 
good initial guess. 
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B. Multi-Dimensional Microwave Sub-wavelength Imaging 


(1) Microwave Sub-wavelength Imaging for Transverse and 
Longitudinal Distributions 
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Fig. 10. Schematic for microwave sub-wavelength imaging using PC lens 
with the objects located along (a) transverse axis (x-axis), (b) longitudinal 
axis (z-axis). 


Figure 10 shows the representative scheme for imaging 
systems using phase conjugating lenses with transversely and 
longitudinally distributed targets. The test objects are PEC 
squares that are illuminated by an incident plane wave. At the 
measurement plane, the phase conjugation lens collects the 
scattered field; phase conjugates this field and propagates it to 
the image plane to form the image. 


For the transverse case, the distance between the source and 
measurement planes is fixed to d, and the distance between 
the image and measurement planes is also fixed to d (in Fig. 
10(a)). However, for the longitudinal case, the PEC squares 
are located along the z-axis, and the image plane is shifting 
along the z-axis in the image region (in Fig. 10(b)). 


(2) SVD-based Basis Matrix Method 


Next, we describe an imaging algorithm based on the 
Singular-Value-Decomposition (SVD)-based basis matrix 
method. As a starter, we use the FEKO 6.2 Suit to simulate the 
case of a single PEC square at a time, which is located in the 
source plane. We derive the scattered field distribution in the 
measurement plane, phase conjugate it, and then let it 
propagate to the image plane. Finally, we extract a part of the 
image distribution, for instance its values along the y = O cut 
to form a column of the basis matrix. The number of the 
columns in this matrix corresponds to the number of possible 
positions for the PEC square in the source plane. We then 


apply a SVD-based basis matrix method to it, and use a 
threshold to delete the singular vectors, whose singular value 
fall below the threshold. 


The step-by-step procedure can be summarized as follows: 


1) Locate one PEC square in the “possible” position at 
one time, and extract a part of the image distribution, for 
instance its values along the y = O cut. 


2) Create the basis matrix by using the extracted values 
obtained in the previous step, and the number of the columns 
in this matrix corresponds to the number of possible positions 
for the PEC square in the source plane. 


3) The basis matrix is written as Bw xy. Hence, we can 
generate an equation as follows: 


(11) 


where Biy x m represents the basis matrix, and N is the number 
of the extracted values and M is the number of possible 
positions for the PEC square in the source plane. Here ex n 
represents the “measured” values, and x,y x 7) represents the 
possible targets. 


En = B ysy ‘X Mx 


4) Equation (11) usually represents an ill-posed problem. 
Hence, we apply the SVD to the basis matrix, 1.e., to derive B 
= UDV, and use a threshold level to reduce the matrix by 
deleting the singular vectors whose singular values fall below 
the threshold. 


5) Take a dot product of the K (K < M) remaining singular 
vectors Uw x x with the “measured” values, as well as the 
basis matrix, to obtain a reduced vector e"x, n and a reduced 
matrix BK xu): 


(12) 


6) Unlike (11), (12) represents a well-conditioned 
problem, and the size of the associated basis matrix B” x m) iS 
significantly reduced compared to that in (11). Hence, we 
directly pseudo-inverse the reduced basis matrix B” xm and 
take its product with the reduced vector e"(x x , to derive the 
final results. 


n n 
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Fig. 11. Raw image (field distribution) results obtained by using the PC 
lens when two objects are located in the transverse plane along the x-axis, and 
their spacing is varied from 24/5 to 94/10, with a 4/10 step. 


To simulate the transverse case, we set the size of the PEC 
square to be 4/20 X 4/20, and the spacing to be 4/4 in both the 
x- and y-directions. The distance between the source and 
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measurement planes is 31. The size of the measurement 
aperture is 304 X 30/1 with 4/20 interval in both x- and 
y-directions. Also the size of the image zone 1s 104 X 104 with 


24/4 interval in both x- and y-directions. The operating 
frequency is 1 GHz, and the plane wave is normally incident. 
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Fig. 13. Imaging results for the 1D case (x-direction) via the SVD-based 
basis matrix method. 
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Fig. 14. Imaging results for the 2D case (x-y plane) via the SVD-based 
basis matrix method. 


Figure 11 shows the raw image (field distribution) results 
obtained by using the PC lens when two objects are located in 
the transverse plane along the x-axis, and their spacing is 
varied from 24/5 to 94/10, with a 4/10 step. From Fig. 11, it is 
obvious that with the increase of the spacing, the two PEC 
squares can be separated more easily. 


Figure 12 shows the imaging models for the 1D and 2D 
cases. The basis matrices are comprised of 13 and 39 columns, 
respectively. The source zone in the x-direction ranges from 
-154/10 to 154/10. Figs. 13 and 14 show the recovered results 
for the 1D and 2D cases; respectively, the stars and the bars 
represent the original and the recovered positions of PEC 
squares; respectively, and the horizontal line at the level of 0.5 
(normalized) is used to distinguish the real objects from the 


false ones. Here, we have set a threshold of 10° for the 
truncation of the singular values. For the 1D case (x-direction), 
we only extract y = 0 cut to generate the basis matrix, and that 
find 13 singular vectors survive. For the 2D case (x-y plane), 
we use three cuts e.g., at y = (-1, 0, 1) x/4/10, and 35 singular 
vectors are survived. From Figs. 13 and 14, we see that the 
basis matrix method applied in conjunction with the SVD is 
capable of accurately recovering the position of the object. 
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Fig. 15. Raw image (field intensity) results obtained by using the PC lens 
when two objects are located along the longitudinal axis (z-axis): region from 
1A to 24 with a 4/10 step. 
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Fig. 16. Image results (field distribution) obtained by using a PC lens when 
two objects are located along the longitudinal axis (z-axis) and the spacing 
between them is varied: (a) raw images (field distribution) before processing; 
(b) results after processing the raw images by using the proposed technique. 


For the longitudinal case, we set the size of the PEC square 
to be 4/20 X 4/20, and let the spacings in the z-direction to 
range from 34/10 to 24. The distance between the source and 
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measurement planes is 31. The size of the measurement 
aperture is 304 X 304 with 4/20 interval in both x- and 
y-directions. And the size of the image zone is 104 X 104 with 
À/4 interval in both x- and y-directions. The operating 
frequency is 1 GHz, and the plane wave is normally incident. 


Figure 15 shows that the longitudinal case is quite different 
from the transverse case. We keep the position of one PEC 
square unchanged, and shift the other PEC square to the 
negative longitudinal axis is from 1/ to 2A with a 4/10 step. It 
shows that the raw image of the two PEC squares merges 
twice with an increase of the spacing when the spacing is 1.4/ 
and 2/, while the image retains its separation in the transverse 
case. Fig. 16 shows the recovered results for the longitudinal 
case—the bars represent the recovered positions of PEC 
squares and we can see that the SVD-based basis matrix 
method also works in the longitudinal case. 


It is evident from the results presented above that the 
SVD-based basis matrix method can enhance the resolution of 
the images generated by a PC lens quite considerably, both in 
the longitudinal and transverse separation cases. 


Ill. LUNEBURG LENS ANTENNAS FOR 
DIRECTION-OF-ARRIVAL ESTIMATION 


A. Direction-of-Arrival Estimation with Normal Luneburg 
Lens Antenna 


The Luneburg lens is designed to focus a plane wave, 
arriving from an arbitrary direction, at a point which is located 
diametrically opposite to that of the incident side. Luneburg 
has shown that the problem of finding €, of the lens can be 
formulated in terms of an integral equation [25], whose 
solution provides us the required material parameters of the 
lens. The radial variation of €, is 


€.=2-(r/RY 


where r is the distance from the center of the lens, and R is its 
radius. 


(13) 


Next, we introduce the principles of sparse reconstruction 
algorithms to analyze the normal Luneburg lens antenna. 
Without loss of generality, we consider a 1D DOA case, in 
which there is only one degree of freedom for the incident 


angle. To realize it, we fix the incident angle g and only vary 0. 


We distribute the sensors on the curved surface of the 
Luneburg lens in a uniform manner and use these sensors to 
measure the electric field at their locations. According to the 
concepts of the sparse signal representation, we generate the 
relationship between the unknown incident angle vector x (M 
x 1) and the measured data vector y (N x 1), where M is the 
number of possible incident angles and N is the number of 
sensors. The matrix A (N x M, N < M) is comprised of training 
data, which is collected or simulated for each possible 
incident angle and is known a priori. Hence, the DOA 
problem with additive white Gaussian can be stated as: 

y=Ax+n (14) 
where n (N x 1) represents the noise vector. 


The /; -norm objective function for (14) can be written as: 


min |x| subjectto |ly — Ax| <p (15) 


where p is a parameter specifying an allowable noise level. 


We use a commercial Finite Difference Time Domain 
(FDTD) [41] code to simulate the antenna geometry. We set 
the interval between the sensors to be 5°, and we arrange 37 
(N = 37) sensors on the surface of the Luneburg lens. The 
number of possible incident angles is set to be 361 (M = 361, 
-90° < 0 < 90°). Since the Luneburg lens is spherically 
symmetric, we only simulate the problem for one incident 
angle, collect the electric field and rotate it by a 0.5° step (A0 
= 0.5°). The incident plane waves are assumed to arrive from 
7 different directions. Their incident angles and normalized 
coefficients are: (62.6°, 0.486, (-6.3°, 1.00e 7 ), (-9.4°, 
0.77e/*), (-44.4°, 0.20e"), (-49.7°, 0.17e%*), (-54.2°, 
0.51e1% and (-71.3°, 0.83e! ). We assume the frequency to 
be 8 GHz though it can be easily changed, and we set the SNR 
to be 30 dB. We apply the algorithm for large-scale sparse 
reconstruction (SPGL1) [42] to the measured data for the 
l -norm sparse reconstruction. The recovered results obtained 
by using the CM and the SPGLI1 algorithms are shown in Fig. 
17. 
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Fig. 17. Field distribution and recovered results by using CM and the 
SPGL1 algorithms. 


Figure 17 shows that the /;-norm sparse reconstruction 
algorithm, 1.e., SPGLI algorithm, has much better resolution 
than the CM. The highest resolution achieved by using 
l—norm sparse reconstruction algorithm is 3.1°, which is 
better than the angular resolution (3.5° @ 8 GHz) predicted on 
the basis of the 3dB-beamwidth concept. Although there is a 
ghost at -75° with a relative small magnitude (around -23 dB), 
considering the SNR is 30 dB, a ghost with such a low level is 
acceptable. 


From Fig. 17, we also see that the proposed method can be 
utilized for wide angle scanning, even to -71.3°, and works for 
ultra-wideband frequency, range of 2 GHz to 12 GHz, 
although we have only shown the 8 GHz case. 


B. Direction-of-Arrival Estimation with Flat-base Luneburg 
Lens Antenna 


(1) Design of Flat-base Luneburg Lens Antenna 


In [29], we have proposed a design, shown in Fig. 18, of a 
lens with a spherical profile, which preserves the salutary 
features of the conventional Luneburg lens but is much more 
convenient to feed in the transmit mode since it has a flat-base 
design. In this section, we consider its advantages in DOA 
estimation by using signal processing techniques. 
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As shown in Fig. 18, the Luneburg lens consists of 11 layers 
and has a diameter of 63.5 mm. The first ten layers from the 
center have a thickness of 3 mm each and the last layer is 1.75 
mm thick. The dielectric parameters of the layers vary 
depending on their distance from the center, in accordance 
with (13). The flat-base is a 6x6 waveguide array with PEC 
walls that are 1 mm thick. The waveguide has a square 
cross-section with a side length of 9 mm. We set an array of 
observers in the measurement plane, 1.e., the truncation plane 
of the waveguides. 


Luneburg Lens 






Waveguide 


Array So 


Measurement 
Plane 


Fig. 18. Flat-base Luneburg lens antenna for DOA estimation. 


It is well-known that the angular resolution of an antenna is 
determined by its half-power beam-width, which is inversely 
proportional to the aperture size of the antenna. We calculate 
the half-power beam-width of the Luneburg lens as A/D by 
assuming the aperture efficiency is unity, where / is the free 
space wavelength and D is the diameter of the Luneburg lens. 
This is slightly larger than 9° for a diameter of D = 63.5 mm at 
30 GHz. 


(2) Correlation Method 


To enhance the resolution of the above Luneburg lens, we 
propose to employ the Correlation Method (CM). In the 
following section we describe the case based on CM to detect 
multiple targets with different incident angles and different 
Radar Cross Section (RCS) levels. We develop the 
measurement database, comprising of the electric fields in the 
measurement plane induced by targets that have different 
combinations of orientations, distances and scattering 
cross-sections. Thus, the measurement database is derived by 
doing a vector addition of the field distributions due to the 
individual targets located at different angles. 

The CM is based on correlating the measured field 
distribution in the measurement plane due to an unknown 
distribution of targets with all the entries in the measurement 
database. The correlation coefficient (p,) for the n-th entry in 
the database is defined as: 


1 
DENE” Cy) -uE y) Y] 


i=l j=l 


Pn 
(16) 


nN 


where the superscript ‘M’ denotes the field distribution in the 
measurement plane induced by the targets; the superscript ‘ D’ 
denotes the known field distribution in the measurement plane, 


which is stored in the database; and u and ø refer to the mean 
and standard deviation of the field distribution; respectively. 
Also, M; and M; are the number of sample points in the 
measurement plane along x- and y-directions, respectively. 
The value of n that maximizes p determines the angular 
locations of the targets, as well as their relative RCSs. 


To demonstrate the application of the method described, we 
now consider the cases of one or two targets, though the 
approach can be extended to detect more than two targets 
located in close proximity. The CM can achieve a very high 
resolution if one can accurately locate the maximum 
correlation coefficient. But, the presence of noise can 
seriously impair the ability to do so if the number of field 
sampling points in the measurement plane is small. Thus, not 
unexpectedly, noise is an important factor that limits the 
resolution. 


We also note that there is a cost involved in terms of higher 
computational complexity, which arises because we need to 
compute the correlation for a large number of combinations. 
For example, if there are two targets within a beam-width and 
they are randomly located in one of the N; possible angles, the 
possible scales of the amplitudes from the two targets is N3, 
the possible initial phase difference between the scattering 
waves from the two targets is N3 then the number of 
combinations which need to considered equals N, xN, xN2xN3. 
which can be rather large. However, the computational burden 
could be reduced by recognizing the fact that that each 
correlation operation is independent and, therefore, the entire 
process is highly parallelizable. Moreover, the scaling 
efficiency in parallel or distributed computations would 
approach unity because the correlation operations with 
different combinations are totally independent and the 
communication between the processes is rarely required. 
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Fig. 19. E., E, and E, field distributions in the measurement plane for 
different incident angles: (a) 9 = 30°; and (b) 0 = 40°, with g = 0°. 


The scattered field from a target which is being interrogated 
by the antenna operating in the receive mode can be treated as 
a plane wave incident up on the antenna from the direction of 
the target. We use a commercial Finite Difference Time 
Domain (FDTD) code to simulate the antenna geometry, 
shown in Fig. 18, to obtain the field distribution in the 
measurement plane. We do this for a y-polarized unit 
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amplitude plane wave for 16 different incident angles, given 
by g = 0° and @ varying from 30° to 45° with an interval of 1°, 
and we set the observation frequency to 30 GHz. Fig. 19 
shows the F,, E,, and E, field distributions at the measurement 
plane for a few selected incident angles. We note that the level 
of E, is always higher than £E, and E, because of the incident 
wave is y-polarized. Moreover, these field distributions vary 
slightly with the gradual change in the angle of incidence. The 
maximum of the field distribution shifts along the negative 
x-direction as 0 varies from 30° to 45°. 


The simulated measured field provides us the positions of 
the samplings at arbitrarily locations in the measurement 
plane, depending upon the mesh used in the simulation. For 
plane waves with different incident angles, the meshes used 
are different and non-uniform. Therefore, we use a linear 
interpolation algorithm to map the field to uniform grids for 
the purpose of generating the database that would be used for 
correlation at a later stage. 
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Fig. 20. Correlation coefficients for (a) E, (b) E, (c) E; and (d) product of 
the three components, for g = 0° and @ varying from 30° to 45°. 
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Fig. 21. Correlation coefficients for a single target located at an angle of 0 
= 37°, g = 0°, with: (a) SNR = 35 dB; (b) SNR = 25 dB; (c) SNR = 15 dB; and 
(d) SNR =5 GB. 


In Fig. 20, the correlation coefficients for each individual 
electric field component, as well as a product of the 
correlation coefficients of all the three components are 
presented. The diagonal elements represent the 
self-correlation coefficients which would be 1 for the ideal 
case when there is no noise. By comparing Fig. 20 (a), (b), (c) 
with (d), it is observed that taking the product of the three 
correlation coefficients further improves the resolution as 
compared to the case when just a single electric field 


component is used. But this comes at the cost of measuring all 
the three components of the electric field, which requires three 
separate sensors in a practical implementation. 


Figure 21 shows the correlation coefficients for a single 
target located at an angle of 0 = 37°, ọ = 0° with different 
noise levels, e.g., SNR = 35 dB, SNR = 25 dB, SNR = 15 dB, 
and SNR = 5 dB. Here, we consider the additive Gaussian 
white noise that has a constant spectral density and a Gaussian 
distribution of amplitude for convenience. We can find that 
the correlation coefficients decrease with an increase of the 
noise level. It means that the results of the correlation method 
are affected by noise. 


We next investigate the case of two targets with very small 
angular separation, 1.e., less than half-power beam-width. The 
two angles were randomly chosen from 30° to 45° in 6 with o 
= (0°, and their combined electric field distributions in the 
measurement plane were derived by using random amplitude 
scaling factor in the range of [0.4, 2.5] and random initial 
phases in the range of [0°, 359°]. We consider 50 random pairs 
of incident plane waves. If there is only one target it would 
correspond to the case when 1“ and 2" angles are the same. 
The total fields measured in the measurement plane can be 
expressed as: 

E,,, ,=ae"E, +aeE, 


total 


| (17) 
= ae” G 4 rete, 


1 


where E; and EE; are the known electrical fields measured in 
the measurement plane when uniform plane waves with unit 
amplitudes are incident on the lens from two different 
directions. Also, (a,, 4) and (a2, Ÿ,) are the amplitude scales 
and initial phases of the fields scattered by the first and second 
randomly chosen scatterers, respectively. 


In the CM approach, the field distributions in the 
measurement plane are then correlated with those in the 
measurement database. The highest value of correlation 
coefficient estimates the angular location of the targets and 
their characteristics. In this paper we consider the step of the 
amplitude scaling factor (a2/a;) as: 1, 1.25, 1.5, 1.75, 2.0, 2.25, 
2.5, 1/1.25, 1/1.5, 1/1.75, 1/2.0, 1/2.25 and 1/2.5; we also 
consider the step of the difference of phase (9>- 9) as: 0°, 9”, 
18°, ..., 351°. We note that, for a pair of incident angles, with 
an amplitude step of 0.25 and phase step of 9°, it is necessary 
to consider a total of 16x16x10x40 different combinations 
that are obtained if we use a combination of 10 cases for the 
amplitude scaling factor (a»/a,), and 40 for the initial phase 
difference (22-81). 


100 


Ratio (%) 
Ratio (%) 
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(18 , 0.25) 
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(9 , 0.5) 


(9 , 0.25) 


(18 , 0.25) (9, 0.5) 
(b) Second angle 

Fig. 22. Accuracy ratios for first (a) and second (b) incident angle with 
different search steps: (i) 69 = 9° and da = 0.25, (ii) 69 = 18° and 6a = 0.25, 
and (iii) 69 = 9° and da = 0.5. 
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Fig. 23. Error distributions and accuracy ratios in the estimation results for 


the first (a) (c) and second (b) (d) incident angle for different SNR with 18 
samplings at each waveguide truncation plane. 
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Fig. 24. Error distributions and accuracy ratios in the estimation results for 
the first (a) (c) and second (b) (d) incident angle for different SNR with 9 
samplings at each waveguide truncation plane. 


Next, we determine the optimal combination of the 
amplitude scale step (6a) and the phase difference step (6@) to 
achieve a high resolution. Three different da and 09 
combinations are considered: (i) 60 = 9° and 6a = 0.25; (ii) 
oO = 18° and da = 0.25; and (iii) dO = 9° and da = 0.5. Fig. 22 
shows the percentage accuracy for these three cases, 
computed by using the 50 random pairs of targets in the 
absence of noise. As expected, the smaller the amplitude scale 
and phase difference steps, the higher the accuracy we can 
achieve. In all cases, the angular error in the detection is less 
than 2°. Also, we note that using (i) results in more accurate 


estimation at the cost of higher computational complexity. 
Hence, we choose 6@ = 9° and da = 0.25 to construct the 
measurement database. 
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Fig. 25. Error distributions and accuracy ratios in the estimation results for 
the first (a) (c) and second (b) (d) incident angle for different SNR with 4 
samplings at each waveguide truncation plane. 
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Fig. 26. Error distributions and accuracy ratios in the estimation results for 
the first (a) (c) and second (b) (d) incident angle for different SNR with 1 
sampling at each waveguide truncation plane. 


Figures 23-26 show the distributions of the target 
detection error for the 50 random locations of a pair of targets 
and their percentage target estimation accuracy, for different 
SNRs and for 18, 9, 4 and 1 field sampling points, 
respectively, in the x- and y-directions for each waveguide in 
the measurement plane, using search steps of 60 = 9° and da 
= 0.25. It is evident that the target estimation accuracy is more 
susceptible to noise when the number of field sampling points 
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is lower. However, since it is difficult to measure the fields at 
a large number of points in a waveguide in practical scenarios, 
we must resort to a compromise between the target detection 
accuracy and the number of sampling points. Also, as the SNR 
decreases to 10 dB, percentage target estimation accuracies 
achieved are 100%, 50%, 40% and 30%, for 18, 9, 4 and 1 
sample, respectively. But, the SNR of a typical radar system 
can be kept well above 20 dB. Therefore, for a SNR above 20 
dB, we are able to achieve percentage target estimation 
accuracies of 100%, 92%, 91% and 65%, respectively, by 
using 18, 9, 4 and 1 sample, respectively. The error in the 
detection in the absence of noise is due to the error introduced 
by interpolation. Also, it is observed from Figs. 23-26 that the 
maximum error, the mean error and the variance of the error 
increase as SNR decreases and the number of field sampling 
point decreases. But these values are within reasonable limits 
for more than 4 field sampling points in the x- and y-directions 
in the measurement plane for each waveguide with SNR 
greater than 20 dB. 


(3) Waveguide Mode Extraction and SVD-based Basis 
Matrix Method 


In section above, a spherical Luneburg lens antenna with a 
planar waveguide array has been proposed. To enhance the 
angular resolution of the lens antenna, we propose a 
waveguide mode extraction technique followed by the use of 
signal processing techniques. We begin by sampling the 
electric field at the cross-section plane of the waveguide array 
and then extracting the weights of the possible propagating 
modes in each waveguide. Next, we apply the CM to obtain an 
initial estimate of the possible DOAs. Following this, we 
employ an SVD-based basis matrix method to achieve a more 
precise estimate of the DOAs. 


Next, we detail the procedure for waveguide mode 
extraction. The cutoff frequencies of the electromagnetic 
waveguide can be complied from [43]: 


1 (mf (ny 
Souto (m,n) = 2 ue =) + z) 


where a and b represent the dimensions of the waveguide 
aperture at the transverse (x-y) plane; € and u are the 
permittivity and permeability of the medium inside the 
waveguide, respectively; and the integers m and n indicate 
different modes. 


(18) 


We set the operating frequency to be f= 30 GHz and choose 
the side dimensions of each square waveguide to be a= b=9 
mm. The cutoff frequencies for TE, 9, TEo;, TE: and TM; 
modes are calculated by using (18) and are found to be 16.67 
GHz, 16.67 GHz, 23.57 GHz, and 23.57 GHz, respectively, as 
shown in Table 1. 

Table 1: Possible propagating modes within different frequency bands. 

Frequency Band 


| gaye Possible Propagating Modes 


(GHz) 


< 16.67 No propagating mode 


16.67 ~ 23.57 TE0 and TEoi 


TE0, TEo1, TE: and 1M.: 
TE 0, TEo1, TE, T™,, and high modes 


2351-3333 
> 33.33 





Once the possible modes are known, the field distribution 
on the face of each waveguide can be expressed as a weighted 
sum of the field distributions for these possible modes as 
follows: 


E=c,E""+c,E°" +¢,E "HE (19) 


where c; (i = 1, 2, 3 and 4) are the weights of each mode and 


ET 10, E7201, ETE11 and E711 are the bases corresponding to 
field distributions for TE; ), TE ),, TE,; and TM,, modes, 
respectively. The field distributions corresponding to the 
possible modes form a complete set of orthogonal bases; 
hence, they can be described any set of field distributions in 
the waveguide under consideration. 


The weights can be related to the sampling data E in matrix 
form as: 


E=A:c (20) 


The next question that needs to be addressed is how many 
field samples are required inside each waveguide to derive an 
accurate DOA estimation. As is known from (19), for the 
structure and frequency range under consideration there are at 
most four degrees of freedom in each waveguide. Viewing 
this issue from the perspective of the rank of the matrix A, the 
maximum rank of A would be 2, 2 and 1, respectively, if we 
only measure the F,-, E,- or E,-component. Hence, we should 
measure both £,- and £,-components to ensure that the matrix 
A is full-rank. 


The next question that needs to be addressed is how many 
field samples are really required inside each waveguide to 
derive an accurate DOA estimation. As is known from (2), for 
the structure and frequency range under consideration there 
are at most four degrees of freedom in each waveguide. 
Viewing this issue from the perspective of the rank of the 
matrix À, if we only measure the E,-, E;- or E,-component, the 
maximum rank of À would be 2, 2 and 1, respectively. Hence, 
we should measure both £,- and £,-components to ensure that 
the matrix A is full-rank. Moreover, if all the sampling points 
are on a Straight line, say along the x- or the y-direction, the 
rank would again reduce to 3; and, hence, a unique solution 
cannot be obtained. To circumvent this problem, we employ 
four samples taken at different x- and y- positions within the 
cross-section of the waveguides. Specifically, we set the 
sample points at (2.25 mm, 2.25 mm), (2.25 mm, 6.75 mm), 
(6.75 mm, 2.25 mm), and (6.75 mm, 6.75 mm) at the 
cross-section plane of each waveguide. 


The step-by-step procedure can be summarized as follows: 


1) Measure the electric field distribution at the planar 
base of the Luneburg lens antenna with a waveguide array for 
different incident angles. Save the observed E,-, and 
F,-components of the field at the sampling points, and 
calculate the weights of the propagating modes in each of the 
36 waveguides for each incident angle. 


2) Create a database matrix using the weights of the 
modes obtained in the previous step. The matrix consists of 
the weights of the modes in all of the waveguides for different 
angles of incidence. 
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3)  Interpolate the database to achieve a finer interval for 
the angle of incidence. The finer database matrix is written as 
Bw». Hence, we can generate an equation as follows: 


=B 


where By x wy represents the database, and N (1.e. 4 x 36 = 144) 
is the number of the weights and M is the number of incident 

angles. Here eyx n represents the “measured” weights, and 

Xm x n represents the possible DOAs. 


e -X (21) 


Nxl NxM M x1 


4) Calculate the weights of the modes from the 
“measured” data obtained from the field distribution due to 
unknown targets. Correlate the “measured” weights, 1.e. e€~vx 1, 
to By x m to determine the most probable target locations. 


5) Extract the columns of the database matrix that 
correspond to the most probable angular regions as a 
sub-database matrix, 1.e., B'iyx1,, where L (L < M) represents 
the possible incident angles in the angular region of the initial 
estimate. Hence, (21) can be rewritten as: 


e,,, =B'vx1-X° 14 (22) 


where the matrix By «1, Corresponds to the most probable 
angular regions, and x°(y x 3) represents the most probable 
regions of the DOA. 


6) Equation (22) usually represents an ill-posed problem. 
Hence, we apply the SVD to the sub-database matrix, 1.e., to 
derive B* = UDV, and use a threshold to reduce the matrix by 
deleting the singular vectors whose singular values fall below 
the threshold. 


7) Take a dot product of the K (K < L) remaining singular 
vectors Uy x x) with the “measured” weights, as well as the 
sub-database matrix, to obtain a reduced weight vector e"(x x 7) 
and a reduced matrix B"(x x 1). 


(23) 


8) Unlike (22), (23) represents a well-conditioned 
problem, and the size of the associated basis matrix B(x x 1) is 
significantly reduced compared to that in (22). Hence, we 
directly pseudo-inverse the reduced basis matrix Br, and 
take its product with the reduced vector e"(x x 7) to derive the 
final results. 


n n S 
e ka =B KkxL'X Lx 


The advantage of carrying out the first step by using the 
CM is to narrow the search region and reduce the 
computational complexity. Once we have determined the 
approximate location of the targets, we can use the second 
step to achieve more accurate estimates of the DOAs. 
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Fig. 27. Comparison of weights rer in two different ways. 


Figure 27 shows a comparison of the weights for the entire 
36 waveguides calculated by using the two methods discussed, 


herein. We see from this figure that the agreement is quite 
good. 


Simulated field 
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Fig. 28. Comparison between simulated E-field and reproduced E-field by 
using waveguide mode extraction. 
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(b) 
Fig. 29. Recovered DOAs by using proposed method. 


Next, we apply the CM to obtain an initial estimate of the 
possible DOAs. Following this we employ the SVD-based 
basis matrix method to achieve a more precise estimate of the 
DOAs. The advantage of carrying out the first step by using 
the CM is to narrow the search region and reduce the 
computational complexity. Once we have determined the 
approximate location of the targets, we can use the second 
step to achieve more accurate estimates of the DOAs. 


Without loss of generality, the structure was excited with 16 
y-polarized plane waves with each a unit amplitude but 16 
different incident angles, corresponding to g = 0° and @ 
varying from 30° to 45° with an interval of 1°. 


Figure 28 shows a comparison between the simulated and 
reproduced E-fields, derived by using the waveguide mode 
extraction technique when 6 = 30°, ọ = 0°. It is evident that 
the waveguide mode extraction technique is able to accurately 
calculate the weights of the propagating modes and reproduce 
the field distribution, with less than 10% difference between 
the two. The 10% error can be attributed to the evanescent 
modes propagating in the waveguide, since the guide is not 
long enough for these modes to decay completely by the time 
they reach the cross-section plane. 


A case of three targets with angular separations as small as 
2° was used to validate the approach. The three incident 
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angles are set to be 31.5”, 34° and 36° in 0, with g = 0°. From 
Fig. 29(a) it is obvious that the inverse matrix method and the 
CM are not able to recover the DOAs. Fig. 29(b) shows that 
the proposed signal processing method can distinguish the 
DOAs with high accuracy even when the targets are located in 
angular proximity of each other. Fig. 29 also shows that the 
recovered DOAs are very close to the real ones; however, the 
recovered scattering cross-sections do not match the real ones, 
although their magnitudes do maintain the same relative order 
w.r.t each other. 


IV. DOA ESTIMATION BY USING APERTURE ANTENNAS 


A. Principles of DOA Estimation by Using Aperture 
Antennas 


Figure 30 shows a representative aperture antenna, namely 
a parabolic reflector. We assume that the sensors receiving the 
signals are distributed in the focal plane of the reflector along 
the y-axis. Our task is to estimate the azimuthal angle of 
arrival by using the SVD-based basis matrix method. 








Sensors \ 


Fig. 30. Geometry of the parabolic reflector with sensors located in its 
focal plane. 


Uniform Linear Array (ULA) of XN sensors receiving I 
signals from the directions @, i = 1,°°:, I. The incident field is 
assumed to be a plane wave originating from a large distance 
away from the antenna. 


The received signal vector can be expressed as follows [33]: 
x(t) = As(t) + n(t) 
x(t) = [x x OT 
A =[a(G,),---,a(@, )] 
s(t) =[5,(1),--45, (OF 
ne) = [nn (OT 


where x(t) represents the signal received at the n-th sensor, 
with n = 1,°*+, N; the matrix A comprises the vector of the 


(24) 


signal waveform; s;(f) is the i-th signal, where i = 1,-°:, J; and 
n,(t) is the additive white Gaussian noise induced in the n-th 
sensor. 


We begin the procedure, proposed herein, by using the 
correlation method as described below. The correlation 
coefficient p between two signals s;(t) and s,(t) is given by: 


E| s OS®) | 
pa =- (25) 
Ells O} JEJ oF | 


where p represents the linear correction (or dependence) 
between the two signals. If p = 1, there is a total positive 
correlation while there is no correlation if p = 0. 





We follow the correlation step, which provides us an initial 
(crude) estimate of the range of the DOA, with the inverse 
matrix method, which is described below. For the noise-free 
case or if the SNR is low, we can use the inverse matrix 
method to find the incident angle of incoming signal as 
follows, by evaluating: 


Ria = Be id (26) 


where Ry is the output of the inverse matrix method. The 
bar-chart representation of the normalized value of R indicates 
the estimates of Z incident angles. In (26), Sy; 1s the received 


signal at the N sensors. B;', is the inverse matrix of the 


IxN 
original base matrix B, which can be derived in advance by 
using / different incident angles. 


The use of our next step is to investigate the SVD-based 
basis matrix method, which is a type of matrix reconstruction 
technique that is useful dealing with ill-conditional matrices, 
which the matrix B may well be for practical situations. We 
note also that many DOA algorithms perform subspace 
decomposition, by using an eigenvalue-eigenvector analysis. 
Ry which is the output of the SVD operator, can be written as 


Ra = Cr, K y (27) 


IxM 


The normalized values show the positions of the incident 


angles, J is the number of the entries of R, C_',, is the inverse 


IxM 
matrix of the new base matrix, which is the truncated version 
of the original base matrix, whose size N has been truncated to 
M by using the SVD, which delete the eigenvectors whose 
eigenvalues fall below a certain threshold, to improve the 
condition number of the new base matrix. We add that the 
Ky; 18 the received signal at the N sensors following the SVD 
projection. 


To summarize, implementation of the proposed DOA 
estimation technique involves the following two steps: 


Step 1: Use correlation method to obtain an initial and crude 
estimate of the DOA. 


Step 2: Utilize inverse matrix method or the SVD-based 
basis matrix method, depending on the SNR level to narrow 
the search range and improve the resolution. 


B. Numerical Results 


For the purpose of the numerical experiments reported 
herein, we set the frequency of operation to be 30 GHz, and 
use a reflector whose focal length (F) is 16.16 mm, and it has 
an aperture diameter (D) of 63.6 mm (6.364). We begin by 
using the HFSS simulator [44] to generate the original base 
matrix for 61 incident angles ranging from —30° to 30°, with a 
1° interval. We assume that there are 67 sensors (0.09544 
separation) on the focal plane along the y-axis; hence N is 67. 
The SVD threshold is set to 10”, and the SNR is assumed to 
be 40 dB. 
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Fig. 31. For the case of simultaneous incidence of two waves coming from 
0° and 3° (a) correlation coefficient; and (b) normalized magnitude. 
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Fig. 32. For the case of simultaneous incidence of two waves coming from 


0° and 3° (a) DOA results obtained by using the entire range; and (b) fine 
search DOA results (normalized) using the proposed algorithm. 
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Fig. 33. For the case of simultaneous incidence of two waves coming from 
0° and 25° (a) correlation coefficient; and (b) fine search DOA results 
(normalized) using the proposed algorithm. 


Figure 31(a) shows the magnitude distribution of the focal 
plane field when two incident waves from 0° and 3° are 
simultaneously incident upon the reflector. From the plot of 
the correlation coefficient shown in Fig. 31(b), we see that we 
can narrow the search range within —2° and +5°. The DOA 
results for this reduced range are shown in Fig. 32(b), and are 
compared with those in Fig. 32(a) when the entire angular 
range is used. We can see that the proposed two-step method 
is demonstrably superior in resolving two incident angles that 
are quite close to each other. Fig. 33(a) plots the correlation 
coefficient. The two angles are 0° and 25°. We note that the 
search region now range from —5° to 30°. Fig. 33(b) illustrates 
the improvement achieved by using the proposed algorithm. 


Note also that the improvement in the resolution obtained 
by using the proposed approach is far superior to that realized 
when using the correlation or other traditional schemes. 


V. DOA ESTIMATION BY USING RANDOMLY DISTRIBUTED 
RADAR ARRAY 


A. Principles of DOA Estimation by Using Randomly 
Distributed Radar Array 


DOA estimation by using the Uniform Circular Arrays 
(UCAs) or the Uniform Linear Arrays (ULAs) is widely used 
in communication and radar fields, and many high angular 
resolution algorithms, e.g., MUSIC, ESPRIT, and etc., have 
been developed for different applications. However, the DOA 
estimation by using the UCA or the ULA has some limitations, 
1.e., strict geometry model requirement, amplitude and phase 
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calibration for the sensors, and unreliability when some 
sensors are not available. In this part, we investigate a 
framework based on the Randomly Distributed Radar Array 
(RDRA) for DOA estimation, and also present a hybrid 
approach which combines the correlation method and the 
SVD-based basis matrix method for achieving much higher 
angular resolution. The hybrid method involves generating a 
data matrix from different incident angles by recording the 
observation for each sensor of the array either by simulation 
or real measurement. Following this, we correlate the 
measurement with generated data matrix, to make an initial 
estimate of the angular regions of possible DOAs. We then 
apply the basis matrix method based on the SVD of the data 
matrix to estimate the angles of incidence accurately. 
Numerical results show that the proposed hybrid approach can 
achieve a high resolution even when the sensors are randomly 
distributed. 
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Fig. 34. (a) Uniform circular array; (b) uniform linear array; and (c) 
randomly distributed radar array. 
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Fig. 35. Mathematic model for DOA estimation by using RDRA. 
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Fig. 36. Schematic for the SVD-based basis matrix method. 


For Fig. 34 to Fig. 36, we show the schematic for UCA; 
ULA; and RDRA, the mathematic model for DOA estimation, 
and the schematic for the SVD-based basis matrix method, 
when using the RDRA. 


B. Numerical Results 


For the purpose of the numerical experiments reported 
herein, we simulate 5 cases with different setting; the 
parameters used are shown in Table 2. For these cases, we 
compare the results within different number of sensors, 
different incident angles, and different type of the array. The 
frequency, SNR, Aperture Size (AS), and 3dB-beamwidth are 
kept unchanged. 


Incident 
angles 
(degree) 


Table 2: Parameters used in simulations for different cases. 
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In Fig. 37, we present the distribution of the sensors (upper 
left), the attenuation and phase delay for each sensor (upper 
middle), the sign-values for the basis matrix (upper right), the 
recovered DOAs by using the correlation method (lower left), 
inverse matrix method (lower middle), and the SVD-based 
basis matrix method (lower right) for 5 different cases. It 
shows that the SVD-based basis matrix method can achieve 
an enhanced angular resolution. Furthermore, the proposed 
framework does not require a strict geometry model and 
works very well even if some sensors are malfunctioned or 
missing. 





Future discussions of this research will focus on how to set 
the phase reference point for the RDRA in practice, as well as 
how to make the echoes from each sensor coherent. 
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VI. CONCLUSION 


In this work, several signal processing techniques; 
namely the CM, the MRPSM, the sparse reconstruction 
method and the SVD-based basis matrix method, have 
been presented and their abilities to resolve at a 
sub-wavelength level for microwave imaging and high 
angular resolution for DOA estimation have been 
examined. We have demonstrated show that by 
combining these signal-processing techniques, the 
performance of various sensing antennas namely the 
phase conjugation lens and Luneburg lens, aperture 
antennas and the randomly distributed radar array can be 
enhanced. We believe that further research in this 
direction would be very productive because of their 
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Fig. 37. Recovered DOAs by using the correlation method, inverse matrix method, and the SVD-based basis matrix method for 5 different cases. 


potential for improving the resolution of sensing 
antennas significantly. 


REFERENCES 


[1] Nieto-Vesperinas, M. and E. Wolf, “Phase conjugation 
and symmetries with wave fields in free space containing 
evanescent components,” J. Opt. Soc. Amer., Vol. 2, No. 
9, pp. 1429-1434, 1985. 

Pendry, J. B. “Negative refraction makes a perfect lens,” 
Phys. Rev. Lett., Vol. 85, No. 18, pp. 3966-3969, 2000. 
Maslovski, S. and S. Tretyakov, “Phase conjugation and 
perfect lensing,’ J. Appl. Phys., Vol. 94, No. 7, pp. 
4241-4243, 2003. 

Rosny, J. de, G. Lerosey and M. Fink, “Theory of 
electromagnetic time-reversal mirrors,” JEEE Trans. 
Antennas Propag., Vol. 58, No. 10, pp. 3139-3149, 2010. 


[2] 
[3] 


[4] 


20 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


[5] Shiroma, G. S., R. Y. Miyamoto, J. D. Roque, J. M. 
Cardenas, and W. A. Shiroma, “A high-directivity 
combined  self-beam/null-steering array for secure 
point-to-point communications,” JEEE Trans. Microw. 
Theory Tech., Vol. 55, No. 5, pp. 838-844, 2007. 

[6] Gaikovich, K. P. “Subsurface near-field scanning 
tomography,” Phys. Rev. Lett., Vol. 98, No. 18, pp. 
183902-183902, 2007. 

[7] Aliferis, I., T. Savelyev, M. J. Yedlin, J.-Y. Dauvignac, A. 
Yarovoy, C. Pichot, and L. Ligfhart, “Comparison of the 
diffraction stack and time-reversal imaging algorithms 
applied to short-range UWB scattering data,” IEEE Int. 
Conf. Ultra-Wideband (ICUWB 2007), Singapore, Sep. 
24-26, 2007. 

[8] Katko, A. R., S. Gu, J. P. Barrett, B. -I. Popa, G. Shvets, 
and S. A. Cummer, “Phase conjugation and negative 
refraction using nonlinear active metamaterials,” Phys. 
Rev. Lett., Vol. 105, No. 12, pp. 123905-123905, 2010. 

[9] Belov, P. A. and Y. Hao, “Sub-wavelength imaging at 
optical frequencies using a transmission device formed 
by a periodic layered metal-dielectric structure operating 
in the canalization regime,” Phys. Rev. B, Vol. 73, No. 11, 
pp. 113110-113110, 2006. 

[10] Eleftheriades, G., and A. Wong, “Holography-inspired 
screens for sub-wavelength focusing in the near field,” 
IEEE Microw. Wireless Compon. Lett., Vol. 18, No. 4, pp. 
236-238, 2008. 

[11] Merlin, R. “Radiationless electromagnetic interference: 
evanescent-field lenses and perfect focusing,” Science, 
Vol. 317, No. 5840, pp. 927—929, 2007. 

[12] Malyuskin, O. and V. Fusco, “Far field subwavelength 
source resolution using phase conjugating lens assisted 
with evanescent-to-propagating spectrum conversion,” 
IEEE Trans. Antennas Propag., Vol. 58, No. 2, pp. 
459—468, 2010. 

[13] Memarian, M. and G. V. Eleftheriades, 
‘“Evanescent-to-propagating wave conversion in 
sub-wavelength metal-strip gratings,’ IEEE Trans. 
Microw. Theory Tech., Vol. 60, No. 12, pp. 3893-3907, 
2012. 

[14] Ge, G.-D., B.-Z. Wang, D. Wang, D. Zhao, and S. Ding, 
“Subwavelength array of planar monopoles with 
complementary split rings based on far-field time 
reversal,” JEEE Trans. Antennas Propag., Vol. 59, No. 1, 
pp. 4345-4350, 2011. 

[15] Katko, A. R., G. Shvets and S. A. Cummer, “Phase 
conjugation metamaterials: particle design and imaging 
experiments,” Journal of Optics, Vol. 14, No. 11, pp. 
114003-1 14003, 2012. 

[16] Park, Y. K. “Subwavelength light focusing and imaging 
via wavefront shaping in complex media,” Progress In 
Electromagnetics Research Symposium (PIERS), 
Guangzhou, China, August 25-28, 2014. 

[17] Sidorenko, P., Y. Shechtman, Y. C. Eldar, O. Cohen, and 
M. Segev, “Sparsity-based sub-wavelength imaging and 
super-resolution in time-resolved and spectroscopic 
instruments,” Progress In Electromagnetics Research 
Symposium (PIERS), Guangzhou, China, August 25-28, 
2014. 

[18] Mittra, R. (Ed.) Computational Electromagnetics - 


Recent Advances and Engineering Applications, Springer, 
New York, chapter 16, pp. 553-574, 2013. 

[19] Gu, X., C. Pelletti, R. Mittra, and Y. Zhang, “Resolution 
enhancement of phase-conjugating lenses by using signal 
processing algorithms,” JEEE Antennas Wireless Propag. 
Lett., Vol. 13, pp. 511-514, 2014. 

[20] Gu, X., C. Pelletti, R. Mittra, and Y. Zhang, “Signal 
processing approach to electromagnetic sub-wavelength 
imaging,” 2013 IEEE Antennas and Propagation Society 
International Symposium (APS/URSI 2013), Orlando, 
Florida, July 7-13, 2013. 

[21]Gu, X., R. Mittra, and Y. Zhang, “Electromagnetic 
sub-wavelength imaging using the basis matrix method in 
conjunction with singular value decomposition (SVD) 
algorithm,” 20/4 IEEE Antennas and Propagation 
Society International Symposium (APS/URSI 2014), 
Memphis, TN, July 6-11, 2014. 

[22] Mittra, R., X. Gu, and Y. Zhang, “Signal processing 
approach to realizing enhanced resolution from imaging 
systems such as lenses, 2014 XXXIth URSI General 
Assembly and Scientific Symposium (URSI/GASS 2014), 
Beijing, China, August 17-23, 2014. 

[23] Balanis, C. A., Modern Antenna Handbook, Wiley, 2008. 

[24] Lafond, O., M. Himdi, H. Merlet, and P. Lebars, “An 
active reconfigurable antenna at 60 GHz based on plate 
inhomogeneous lens and feeders,” IEEE Trans. Antennas 
Propag., Vol. 61, No. 4, pp. 1672—1678, 2013. 

[25]Luneburg, R. K. Mathematical Theory of Optics, 
University of California Press, 1964. 

[26]James, G., A. Parfitt, J. Kot, and P. Hall, “A case for the 
Luneburg lens as the antenna element for the 
square-kilometre array radio telescope,” Radio Science 
Bulletin, No. 293, pp. 32-37, June, 2000. 

[27]Hua, C., X. Wu, N. Yang, and W. Wu, “Air-filled 
parallel-plate cylindrical modified Luneberg lens antenna 
for multiple-beam scanning at millimeter-wave 
frequencies,” IEEE Trans. Microw. Theory Tech., Vol. 
61, No. 1, pp. 436—443, 2013. 

[28] Liang, M., X. Yu, S.-G., Rafael, W.-R. Ng, M. E. Gehm, 
H. Xin, “Direction of arrival estimation using Luneburg 
lens,” IEEE International Microwave Symposium (IMS) 
Digest (MTT), pp.1-3, June 17-22, 2012. 

[29]Jain S. and R. Mittra, “Flat-base broadband multibeam 
Luneburg lens for wide angle scan’, 20/4 IEEE Antennas 
and Propagation Society International Symposium 
(APS/URSI 2014), Memphis, TN, July 6-11, 2014. 

[30] Gu, X., S. Jain, R. Mittra, and Y. Zhang, “Enhancement 
of angular resolution of a flat-base Luneburg lens antenna 
by using correlation method,” Progress In 
Electromagnetics Research M, Vol. 37, pp. 203-211, 
2014. 

[31] Kasparis, C. “Improved resolution DoA estimation 
through shrunk projections on the signal sub-space,” 
Signal Processing, Vol. 92, pp. 1673—1678, 2012. 

[32] Liao, B.,. Z. Zhang, and S. Chan, “DOA estimation and 
tracking of ulas with mutual coupling,” IEEE 


Transactions on Aerospace and Electronic Systems, Vol. 
48, pp. 891—905, 2012. 


21 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


[33] Gu, X. and Y. Zhang, “Resolution threshold analysis of 
MUSIC algorithm in radar range imaging,” Progress In 
Electromagnetics Research B, Vol. 31, pp.297-321, 
2011. 


[34] Rouquette, S. and M. Najim, “Estimation of frequencies 
and damping factors by two-dimensional ESPRIT type 
methods,” JEEE Trans. Signal Processing, Vol. 49, No. 1, 
pp.237-245, 2001. 

[35] Mittra, R., C. Pelletti, N. L. Tsitsas, and G. Bianconi, “A 
new technique for efficient and accurate analysis of FSSs, 
EBGs and metamaterials,” Microw. Opt. Techn. Lett., 
Vol. 54, No. 4, pp.1108-1116, 2011. 

[36] Pelletti, C., G. Bianconi, R. Mittra, A. Monorchio, and K. 
Panayappan, “Numerically efficient method-of-moments 
formulation valid over a wide frequency band including 
very low frequencies,” JET Microw. Antennas Propag., 
Vol. 6, No. 1, pp. 46-51, 2012. 

[37] FEKO Suite 6.2 [Online]. Available: www.feko.info. 

[38] Donoho, D. L. “Compressed sensing,” IEEE Trans. Inf. 
Theory, Vol. 52, No. 4, pp. 1289-1306, 2006. 

[39] Balanis, C. A. Antenna Theory - Analysis and Design, 
Second edition, John Wiley & Sons, 1982. 

[40] Mohimani, H., M. Babaie-Zadeh, and C. Jutten, “A fast 
approach for overcomplete sparse decomposition based 
on smoothed LO norm,” JEEE Trans. Signal Process., 
Vol. 57, No. 1, pp. 289-301, 2009. 

[41] Yu, W., X. Yang, Y. Liu, R. Mittra, and A. Muto, 
Advanced FDTD Methods: Parallelization, Acceleration, 
and Engineering Applications, Artech House, Norwood, 
MA, USA, Mar. 2011. 

[42] Berg, E.V.D. and M.P. Friedlander, “Sparse optimization 
with least-squares constraints,” SIAM J. OPTIM., Vol. 21, 
No. 4, pp. 1201-1229, 2011. 

[43] Guru B. and H. Hiziroglu, Electromagnetic Field Theory 
Fundamentals, Second edition, Cambridge University 
Press, 2004. 


[44] http://www.ansys.com/ 


APPENDIX 


Phase Conjugating 


MRPSM Minimum Residual Power Search 
Method 


SD | Singular-Vale-Decomposton 
(DOA | Direcion-oFArival 
PTRM | Time Reversal Mitros 
Owe [Drewen 


ESPRIT Estimation of Signal Parameters via 
Rotation Invariance Techniques 


Perfect Electric Conductor 


SVD 
DOA 
TRM 
UWB 
EC 





SNR 
MRPC 


Signal-to-Noise Ratio 


Minimum Residual Power Criterion 


SMBFs Sinusoidal Macro Basis Functions 
SRBFs Sinusoidal Rooftop Basis Functions 
Compressive Sensing 


Least Square Method 


SLO Smoothed lo algorithm 


Finite Difference Time Domain 


Radar Cross Section 
Aperture Size (AS) 





Xiang Gu was born in 1982. He received the 
B.S. degree in the Department of Electrical 
Engineering from Beijing Institute of 
Technology (BIT), Beijing, China, in 2004, 
and the Ph.D. degree from the Graduate 
School of the Chinese Academy of Sciences, 
Beijing, China, in 2009. He was a research 
assistant with the Electromagnetic Communication 
Laboratory (EMC), Department of Electrical Engineering of 
the Pennsylvania State University (PSU). He is an associate 
professor with the Department of Microwave Remote Sensing 
and Information Technology, Center for Space Science and 
Applied Research (CSSAR), Chinese Academy of Sciences. 
His research interests include microwave imaging, 
computational electromagnetics, remote sensing and radar 
signal processing. 





22 


